function Master_mixed_FEdemoshrt_cHD_wtp_varagin_j_bs(subsample,inc,set_seed,matlab_seed,get_nb_draws, j, index_bs)


tic

global state_i week_i J T S MT MTvector K WB...
price_num rebate_estar_num kwh_num elec_num estar_num ...
price_denom rebate_estar_denom rebate_denom kwh_denom elec_denom estar_denom ...
bchoice bchoice_s bchoice_ij id_i id_j id_ij id_m id_n  demo pidnest choicenest_ij brandweek_denom brandweek_num prod_denom prod_num...
type_id2_3_demo1_num type_id2_3_demo1_denom type_id3_demo1_num type_id3_demo1_denom AV_demo1_num AV_demo1_denom  ice_sc_demo1_num ice_sc_demo1_denom...
type_id2_3_demo2_num type_id2_3_demo2_denom type_id3_demo2_num type_id3_demo2_denom AV_demo2_num AV_demo2_denom  ice_sc_demo2_num ice_sc_demo2_denom...
type_id2_3_demo3_num type_id2_3_demo3_denom type_id3_demo3_num type_id3_demo3_denom AV_demo3_num AV_demo3_denom  ice_sc_demo3_num ice_sc_demo3_denom...
type_id2_3_demo4_num type_id2_3_demo4_denom type_id3_demo4_num type_id3_demo4_denom AV_demo4_num AV_demo4_denom  ice_sc_demo4_num ice_sc_demo4_denom...
type_id2_3_demo5_num type_id2_3_demo5_denom type_id3_demo5_num type_id3_demo5_denom AV_demo5_num AV_demo5_denom  ice_sc_demo5_num ice_sc_demo5_denom...
type_id2_3_demo6_num type_id2_3_demo6_denom type_id3_demo6_num type_id3_demo6_denom AV_demo6_num AV_demo6_denom  ice_sc_demo6_num ice_sc_demo6_denom;

global id; 

%Parameters for the mixed logit:
global N_id_denom;
global seq1_Ndraws seq2_Ndraws;
global nb_draws;
global max_size;
global id_dense_i id_dense_j id_dense_ij id_dense_m id_dense_n;
global bchoice_R;

global LLC_5_18;


nb_draws=get_nb_draws;
burn_draws=15;

%Create Diary
date_now = clock;
date_now = strcat(num2str(date_now(1)),'_',num2str(date_now(2)),'_', num2str(date_now(3)), num2str(date_now(4)), num2str(date_now(5)));


%Create Output: See below and make changes
results_file1=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Results_Estimation\','mixed_logit_FEdemoshrt_cHD_wtp_',num2str(subsample),'_',num2str(inc),'_',num2str(set_seed),'_matlab_seed_',num2str(matlab_seed),'_nb_draws_',num2str(nb_draws),'bs_',num2str(j),'.mat');

%Read Files
	tmp_file1=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    id=load(tmp_file1); 
    id=id(index_bs,:);
    id(:,1:7)=[];
    choicest_size=sum(id,2);
    %id_discard=find(choicest_size>200 | choicest_size<50);
    id_discard=find(choicest_size>400);
    id(id_discard,:)=[];  
    choicest_size(id_discard)=[]; 
    id=id';

    MT=size(id,2);
    MTvector=ones(MT,1);
    J=size(id,1);

    %Using information on the size bought
    %we restrict the consideration set
    tmp_file_sizebought=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation/sizebought_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    sizebought_tmp=load(tmp_file_sizebought); 
    sizebought_tmp=sizebought_tmp(index_bs,:);
    sizebought_tmp(:,1)=[];
    sizebought=repmat(sizebought_tmp,1,J);
    sizebought=sizebought';
    
    tmp_file_sizeoffer=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation/sizeoffer_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    sizeoffer_tmp=load(tmp_file_sizeoffer); 
    sizeoffer_tmp(:,1)=[];
    sizeoffer=repmat(sizeoffer_tmp,1,MT);
   
    %The consideration set is defined as all
    %refrigerators withing 2.5 cu. ft. of the size 
    %bought
    sizediff_tmp1=abs(sizebought-sizeoffer);
    sizediff_tmp2=find(sizediff_tmp1<2.5);
    sizediff_tmp3=zeros(J,MT);
    sizediff_tmp3(sizediff_tmp2)=1;
    id_size=id.*sizediff_tmp3;
    choicest_size_cHD=sum(id_size,1);
    clear sizebought_tmp sizebought sizeoffer_tmp sizeoffer sizediff_tmp1 sizediff_tmp2 sizediff_tmp3;
    
    %id_size=id_size';
    id_R=repmat(id_size,1,nb_draws);
    [id_i,id_j,id_s]=find(id_R);
    id_ij=find(id_R);
    [id_m,id_n] = size(id_R);
    
    MTR=size(id_R,2);
    %clear id_size id_R;
    
    max_size=max(choicest_size_cHD);
    for i=1:MT
     id_dense(i,:)=[ones(1,choicest_size_cHD(i)),zeros(1,max_size-choicest_size_cHD(i))];
    end
    id_dense=id_dense';
    id_dense_R=repmat(id_dense,1,nb_draws);
    
    [id_dense_i,id_dense_j,id_dense_s]=find(id_dense_R);
    id_dense_ij=find(id_dense_R);
    [id_dense_m,id_dense_n] = size(id_dense_R);
    clear id_dense id_dense_R;
    
 	%bchoice=load('H:\Research\sears\estar_data\refrigerators\bchoice_046_2008_2012_v11022017_struct_52171_4.csv'); 
    tmp_file2=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\bchoice_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
 	bchoice=load(tmp_file2);
 	bchoice=bchoice(index_bs,:); 
    bchoice(:,1)=[];
    bchoice(id_discard,:)=[]; 
    bchoice=bchoice';
    bchoice_R=repmat(bchoice,1,nb_draws);
    [bchoice_i,bchoice_j,bchoice_s]=find(bchoice_R);
    bchoice_ij=find(bchoice_R);
    [bchoice_m,bchoice_n] = size(bchoice_R);
    clear bchoice;    
    
    tmp_file3=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\estar_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    estar=load(tmp_file3); 
    estar=estar(index_bs,:);
    estar(:,1:7)=[];
    estar(id_discard,:)=[]; 
    estar=estar';
    estar_R=repmat(estar,1,nb_draws);
    estar_num=estar_R(bchoice_ij);
    estar_denom=estar_R(id_ij);
    clear estar estar_R;
    
    tmp_file4=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\rebate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    rebate=load(tmp_file4); 
    rebate=rebate(index_bs,:); 
    rebate(:,1:7)=[];
    rebate=rebate/1000;
    rebate_J=repmat(rebate,1,J);
    rebate_J(id_discard,:)=[]; 
    rebate_J=rebate_J';
    rebate_J_R=repmat(rebate_J,1,nb_draws);  
    rebate_denom=rebate_J_R(id_ij);
    rebate_estar_num=rebate_J_R(bchoice_ij).*estar_num;
    rebate_estar_denom=rebate_J_R(id_ij).*estar_denom;
    clear rebate rebate_J rebate_J_R;
    
    tmp_file5=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\tax_rate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tax=load(tmp_file5); 
    tax=tax(index_bs,:);
    tax(:,1:7)=[];
    tax(id_discard,:)=[]; 
    taxp=1+repmat(tax',J,1);
    
    tmp_file6=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\price_mean_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    price_net=load(tmp_file6); 
    price_net=price_net(index_bs,:);
    price_net(:,1:7)=[];
    price_net(id_discard,:)=[]; 
    price_net=price_net';
    price=taxp.*price_net;
    price=price/1000;
    price_R=repmat(price,1,nb_draws); 
    price_num=price_R(bchoice_ij);
    price_denom=price_R(id_ij);
    clear price_net taxp tax price_R;
    
    tmp_file7=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\kwh_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    kwh=load(tmp_file7); 
    kwh=kwh(index_bs,:);
    kwh(:,1:7)=[];
    kwh(id_discard,:)=[]; 
    kwh=kwh';
    kwh=kwh/1000;
    kwh_R=repmat(kwh,1,nb_draws); 
    kwh_num=kwh_R(bchoice_ij);
    kwh_denom=kwh_R(id_ij);
    clear kwh kwh_R;
    
    tmp_file8=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\elec_county_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    elec=load(tmp_file8); 
    elec=elec(index_bs,:);
    elec(:,1:7)=[];
    elec_J=repmat(elec,1,J);
    elec_J(id_discard,:)=[]; 
    elec_J=elec_J';
    elec_J_R=repmat(elec_J,1,nb_draws); 
    elec_num=elec_J_R(bchoice_ij).*kwh_num;
    elec_denom=elec_J_R(id_ij).*kwh_denom;
    clear elec_J elec_J_R;
    
    tmp_file9=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\demo_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    demo=load(tmp_file9); 
    demo=demo(index_bs,:);

        %hd_id income2 income3 education2 education3 fam_size age political2 political3
    demo(:,1)=[];
    demo(:,1:2)=[];
    demo(:,3)=demo(:,3)/10;
    demo(:,4)=demo(:,4)/100;
    demo1_J=repmat(demo(:,1),1,J); %education2
    demo2_J=repmat(demo(:,2),1,J); %education3
    demo3_J=repmat(demo(:,3),1,J); %fam_size
    demo4_J=repmat(demo(:,4),1,J); %age
    demo5_J=repmat(demo(:,5),1,J); %political2
    demo6_J=repmat(demo(:,6),1,J); %political3
   
%Create Interactions Attributes-Demographics     
	tmp_file10=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\type_id_2_3_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
	type_id2_3=load(tmp_file10); 
	type_id2_3=type_id2_3(index_bs,:);
    type_id2_3(:,1:7)=[];
   
    type_id2_3_demo1=type_id2_3.*demo1_J;
    type_id2_3_demo1(id_discard,:)=[]; 
    type_id2_3_demo1=type_id2_3_demo1';
    type_id2_3_demo1_R=repmat(type_id2_3_demo1,1,nb_draws);
    type_id2_3_demo1_num=type_id2_3_demo1_R(bchoice_ij);
    type_id2_3_demo1_denom=type_id2_3_demo1_R(id_ij);
    clear type_id2_3_demo1 type_id2_3_demo1_R;
    
    type_id2_3_demo2=type_id2_3.*demo2_J;
    type_id2_3_demo2(id_discard,:)=[];
    type_id2_3_demo2=type_id2_3_demo2';
    type_id2_3_demo2_R=repmat(type_id2_3_demo2,1,nb_draws);
    type_id2_3_demo2_num=type_id2_3_demo2_R(bchoice_ij);
    type_id2_3_demo2_denom=type_id2_3_demo2_R(id_ij);
    clear type_id2_3_demo2 type_id2_3_demo2_R;
    
    type_id2_3_demo3=type_id2_3.*demo3_J;
    type_id2_3_demo3(id_discard,:)=[];
    type_id2_3_demo3=type_id2_3_demo3';
    type_id2_3_demo3_R=repmat(type_id2_3_demo3,1,nb_draws);
    type_id2_3_demo3_num=type_id2_3_demo3_R(bchoice_ij);
    type_id2_3_demo3_denom=type_id2_3_demo3_R(id_ij);
    clear type_id2_3_demo3 type_id2_3_demo3_R;
    
    type_id2_3_demo4=type_id2_3.*demo4_J;
    type_id2_3_demo4(id_discard,:)=[];
    type_id2_3_demo4=type_id2_3_demo4';
	type_id2_3_demo4_R=repmat(type_id2_3_demo4,1,nb_draws);   
    type_id2_3_demo4_num=type_id2_3_demo4_R(bchoice_ij);
    type_id2_3_demo4_denom=type_id2_3_demo4_R(id_ij);
    clear type_id2_3_demo4 type_id2_3_demo4_R;
    
    type_id2_3_demo5=type_id2_3.*demo5_J;
    type_id2_3_demo5(id_discard,:)=[];
    type_id2_3_demo5=type_id2_3_demo5';
   	type_id2_3_demo5_R=repmat(type_id2_3_demo5,1,nb_draws);
    type_id2_3_demo5_num=type_id2_3_demo5_R(bchoice_ij);
    type_id2_3_demo5_denom=type_id2_3_demo5_R(id_ij);
    clear type_id2_3_demo5 type_id2_3_demo5_R;
    
    type_id2_3_demo6=type_id2_3.*demo6_J;
    type_id2_3_demo6(id_discard,:)=[];
    type_id2_3_demo6=type_id2_3_demo6';
    type_id2_3_demo6_R=repmat(type_id2_3_demo6,1,nb_draws);
    type_id2_3_demo6_num=type_id2_3_demo6_R(bchoice_ij);
    type_id2_3_demo6_denom=type_id2_3_demo6_R(id_ij);
    clear type_id2_3_demo6 type_id2_3_demo6_R;
   
   
    tmp_file11=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\AV_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    AV=load(tmp_file11);
    AV=AV(index_bs,:); 
    AV(:,1:7)=[];
    AV=AV/100;
    
    AV_demo1=AV.*demo1_J;
    AV_demo1(id_discard,:)=[];
    AV_demo1=AV_demo1';
    AV_demo1_R=repmat(AV_demo1,1,nb_draws);
    AV_demo1_num=AV_demo1_R(bchoice_ij);
    AV_demo1_denom=AV_demo1_R(id_ij);
    clear AV_demo1 AV_demo1_R;
    
    AV_demo2=AV.*demo2_J;
    AV_demo2(id_discard,:)=[];
    AV_demo2=AV_demo2';
    AV_demo2_R=repmat(AV_demo2,1,nb_draws);
    AV_demo2_num=AV_demo2_R(bchoice_ij);
    AV_demo2_denom=AV_demo2_R(id_ij);
    clear AV_demo2 AV_demo2_R;
    
    AV_demo3=AV.*demo3_J;
    AV_demo3(id_discard,:)=[];
    AV_demo3=AV_demo3';
    AV_demo3_R=repmat(AV_demo3,1,nb_draws);
    AV_demo3_num=AV_demo3_R(bchoice_ij);
    AV_demo3_denom=AV_demo3_R(id_ij);
    clear AV_demo3 AV_demo3_R;
    
    AV_demo4=AV.*demo4_J;
    AV_demo4(id_discard,:)=[];
    AV_demo4=AV_demo4';
    AV_demo4_R=repmat(AV_demo4,1,nb_draws);
    AV_demo4_num=AV_demo4_R(bchoice_ij);
    AV_demo4_denom=AV_demo4_R(id_ij);
    clear AV_demo4 AV_demo4_R;
    
    AV_demo5=AV.*demo5_J;
    AV_demo5(id_discard,:)=[];
    AV_demo5=AV_demo5';
    AV_demo5_R=repmat(AV_demo5,1,nb_draws);
    AV_demo5_num=AV_demo5_R(bchoice_ij);
    AV_demo5_denom=AV_demo5_R(id_ij);
    clear AV_demo5 AV_demo5_R;
    
    AV_demo6=AV.*demo6_J;
    AV_demo6(id_discard,:)=[];
    AV_demo6=AV_demo6';
    AV_demo6_R=repmat(AV_demo6,1,nb_draws);
    AV_demo6_num=AV_demo6_R(bchoice_ij);
    AV_demo6_denom=AV_demo6_R(id_ij);
    clear AV_demo6  AV_demo6_R;
   
    clear AV;
    
    tmp_file12=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\ice_sc_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    ice_sc=load(tmp_file12); 
    ice_sc=ice_sc(index_bs,:); 
    ice_sc(:,1:7)=[];
   
    ice_sc_demo1=ice_sc.*demo1_J;
    ice_sc_demo1(id_discard,:)=[];
    ice_sc_demo1=ice_sc_demo1';
    ice_sc_demo1_R=repmat(ice_sc_demo1,1,nb_draws);
    ice_sc_demo1_num=ice_sc_demo1_R(bchoice_ij);
    ice_sc_demo1_denom=ice_sc_demo1_R(id_ij);
    clear ice_sc_demo1 ice_sc_demo1_R;
    
    ice_sc_demo2=ice_sc.*demo2_J;
    ice_sc_demo2(id_discard,:)=[];
    ice_sc_demo2=ice_sc_demo2';
    ice_sc_demo2_R=repmat(ice_sc_demo2,1,nb_draws);
    ice_sc_demo2_num=ice_sc_demo2_R(bchoice_ij);
    ice_sc_demo2_denom=ice_sc_demo2_R(id_ij);
    clear ice_sc_demo2 ice_sc_demo2_R;
    
    ice_sc_demo3=ice_sc.*demo3_J;
    ice_sc_demo3(id_discard,:)=[];
    ice_sc_demo3=ice_sc_demo3';
    ice_sc_demo3_R=repmat(ice_sc_demo3,1,nb_draws);
    ice_sc_demo3_num=ice_sc_demo3_R(bchoice_ij);
    ice_sc_demo3_denom=ice_sc_demo3_R(id_ij);
    clear ice_sc_demo3 ice_sc_demo3_R;
    
    ice_sc_demo4=ice_sc.*demo4_J;
    ice_sc_demo4(id_discard,:)=[];
    ice_sc_demo4=ice_sc_demo4';
    ice_sc_demo4_R=repmat(ice_sc_demo4,1,nb_draws);
    ice_sc_demo4_num=ice_sc_demo4_R(bchoice_ij);
    ice_sc_demo4_denom=ice_sc_demo4_R(id_ij);
    clear ice_sc_demo4 ice_sc_demo4_R;
    
    ice_sc_demo5=ice_sc.*demo5_J;
    ice_sc_demo5(id_discard,:)=[];
    ice_sc_demo5=ice_sc_demo5';
    ice_sc_demo5_R=repmat(ice_sc_demo5,1,nb_draws);
    ice_sc_demo5_num=ice_sc_demo5_R(bchoice_ij);
    ice_sc_demo5_denom=ice_sc_demo5_R(id_ij);
    clear ice_sc_demo5 ice_sc_demo5_R;
    
    ice_sc_demo6=ice_sc.*demo6_J;
    ice_sc_demo6(id_discard,:)=[];
    ice_sc_demo6=ice_sc_demo6';
    ice_sc_demo6_R=repmat(ice_sc_demo6,1,nb_draws);
    ice_sc_demo6_num=ice_sc_demo6_R(bchoice_ij);
    ice_sc_demo6_denom=ice_sc_demo6_R(id_ij);
    clear ice_sc_demo6 ice_sc_demo6_R;
   
    clear ice_sc;
      
    prod=repmat([1:J],MTR,1)';
    prod_denom=prod(id_ij);
    prod_num=prod(bchoice_ij);
    
 clear data;   


 
%Initialize Parameters

result_homFE=strcat('Z:\FAC\HEC\DEEP\shoude\default\D2c\eegap\EEgap_data_code_heter_SM\Results_Estimation\hom_FEdemoshrt_cHD_wtp_LLC_5_18_',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_07092023_bs_',num2str(j),'.mat');
eval(['load ' result_homFE]);
result_homFE0=theta_est1;

alpha=result_homFE0(1:J-1,1);
eta=result_homFE0(J)  
tau=result_homFE0(J+1)
pi=result_homFE0(J+2)
theta=result_homFE0(J+3)
beta=result_homFE0(J+4:J+15)  


    %Product FE
        %alpha=ones(1,J-1);
   %Price Elasticity
        %eta=-0.00524720859088;
    %Estar    
        %tau= 0.186990388021096 ;
    %Rebate
        %pi=0.0002;
    %Electricity   
        %theta=-0.01652660097783;
    %Demo 
    	%beta=0.01*ones(1,12);
    %Cov. matrix mixed logit
    	omega=[0.5,01,0.5];
        
  %Homogeneous With Product FEs: ChoiceProb_hom_FEnocol_vs
   %FE
   param_b=[alpha; eta; tau; pi; theta; beta; omega'];
   %Behavioral param
   %param_b(J,1)=eta;
   %param_b(J+1,1)=tau;
   %param_b(J+2,1)=pi;
   %param_b(J+3,1)=theta;
   %AttributesXdemo: 6-type_id2_3, 6-AV
   %param_b(J+4:J+15,1)=beta;
   %param_b(J+16:J+18,1)=omega;
  
   Nparam=J+18;
   
  rho_5=1/(1+0.05);
  LLC_5_18=rho_5*(1-rho_5^18)/(1-rho_5);
toc 

%check_b=ChoiceProb_mixed_FEdemoshrt_grID_wtp_tomlab(param_b);
%results_0=load('\\c3\rdat\SHoude\Research\Sears\estar_results\struct_est\mixed_logit_FEdemoshrt_grID_wtp_11000_16_1_matlab_seed_1_nb_draws_50_bup01222018');
%theta_est0=results_0.theta_est1;
%check_0=ChoiceProb_mixed_FEdemoshrt_grID_wtp_tomlab(theta_est0);

%%Generate Halton sequence

 [seq1_Hdraws,seq2_Hdraws,seq1_Ndraws_MTblocks,seq2_Ndraws_MTblocks]=halton_2normal(burn_draws,nb_draws,MT);
 % The structure of the Halton sequence is in block of observations.
 % [block i=1 [1...R]; block i=2 [1...R]; ...; block i=MT [1...R]];
 % Re-order so that we have the following block structure:
 % [block r=1 [1...MT]; block r=2 [1...MT]; ...; block r=R [1...MT]];
 

  seq1_tmp=reshape(seq1_Ndraws_MTblocks,nb_draws,MT);
  seq2_tmp=reshape(seq2_Ndraws_MTblocks,nb_draws,MT);
  
  seq1_Ndraws=seq1_tmp(:);
  seq2_Ndraws=seq2_tmp(:);
 
% %Using Tomlab-Knitro
% tic
% 	        Name = 'ML Knitro';
% 			
% 			x_L=-1000*ones(Nparam,1);
% 			x_U=1000*ones(Nparam,1);
% 			fLowBnd = -inf;    % Lower bound on function.
% 			tol_outer = 1.e-12;
% 			hess_pat=ones(Nparam);
% 		
% 			%Prob = conAssign('ChoiceProb_mixed_FEdemoshrt_tomlab','ML_grad_hom_FEdemoshrt_tomlab',[],hess_pat, x_L, x_U, Name, param_b,[], fLowBnd, [], [], [], [], [], [], [], [], []);
% 			Prob = conAssign('ChoiceProb_mixed_FEdemoshrt_grID_wtp_tomlab',[],[],hess_pat, x_L, x_U, Name, param_b,[], fLowBnd, [], [], [], [], [], [], [], [], []);
% 			
% 			  tic
%                 [lle]=ChoiceProb_mixed_FEdemoshrt_grID_wtp_tomlab(param_b);
%               toc 
%               
% 			    Prob.PriLevOpt = 3;
% 			    Prob.KNITRO.PrintFile = ''; 
% 			    Prob.KNITRO.options.MAXIT = 500;    % Setting maximum number of iterations     
% 			    Prob.KNITRO.options.OPTTOL = tol_outer ;    % set outer-loop tolerance
% 			      
%                 Prob.KNITRO.options.PRESOLVE = 0; 
% 			    Prob.KNITRO.options.ALG = 1;
% 			    Prob.KNITRO.options.HESSOPT = 2;    % to be used along with first-order derivatives          
% 			    Knitro_Result = tomRun('knitro', Prob, 1);
% 			    theta_est1 = Knitro_Result.x_k
% 			    fvalue_Knitro = Knitro_Result.f_k;
% 			    INFO_Knitro = Knitro_Result.Inform;
%                 Hessian_Knitro=Knitro_Result.H_k;
%                 gradient_Knitro=Knitro_Result.g_k;
% toc_knitro=toc
%eval(['save ' results_file1 ' J Knitro_Result theta_est1 fvalue_Knitro INFO_Knitro toc_knitro Hessian_Knitro gradient_Knitro']);  
 
options_unc4=optimoptions(@fminunc,'Algorithm','trust-region','GradObj','on','DerivativeCheck','off','Display','iter','TolFun',1e-14,'TolX',1e-12,'MaxFunEvals',10000,'MaxIter',2500);

%Using stantard fminunc
tic
  %[theta_est0,fval0,exitflag0,output0]=fminsearch(@(param) lle2opt(param),param_b,options_search);
  [theta_est1,fval1,exitflag1,output1,grad1,hessian1]=fminunc(@(param) ChoiceProb_mixed_FEdemoshrt_grID_wtp_tomlab(param_b),options_unc4);
toc_fminunc=toc

eval(['save ' results_file1 ' theta_est1 fval1 exitflag1 output1 grad1 hessian1 toc_fminunc']);  

     
end 






      





